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We simulate two dynamical, mass degenerate light quarks on 16 3 x 32 lattices with a spatial 
extent of 2.4 fm using the Chirally Improved Dirac operator. The simulation method, the 
implementation of the action and signals of equilibration are discussed in detail. Based on 
the eigenvalues of the Dirac operator we discuss some qualitative features of our approach. 
Results for ground state masses of pseudoscalar and vector mesons as well as for the nucleon 
and delta baryons are presented. 
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T 1 ■ I. INTRODUCTION 

. Lattice Dirac operators that obey the so-called Ginsparg- Wilson (GW) relation [lj implement a 

lattice version of the chiral symmetry transformations [2|. Presently only one explicit formulation 
of lattice fermions, the overlap Dirac operator [1,0], is GW exact in that sense. There are, however, 
several formulations approaching GW exactness in various ways. Among them is the domain-wall 
formulation 0, @] , which approaches the overlap operator in the limit of infinite extent of an artificial 
5-th dimension. Another one is the so-called perfect Dirac operator 0, 0], which if constructed 
explicitly, would obey the GW condition, and which has been approximated by a parameterized 
fixed-point form. Here we discuss a simulation with the so-called Chirally Improved (CI) Dirac 
operator 0, which is also a parameterization of a Dirac operator obeying the GW relation 
■ approximately. 

J> , Advantages of GW exact fermions are that there is no additive mass renormalization and thus no 

spurious zero modes at non-zero quark masses. Operators are protected by chiral symmetry which 
is convenient for the determination of certain matrix elements. Technically the GW exact overlap 
operator involves taking the square root of a simpler kernel operator (e.g., the Wilson operator), 
which is computationally roughly two orders of magnitude more expensive than simulations with 
non-GW-operators. This is not only due to the technical implementation (through, e.g., polynomial 
series or rational functions) but also due to tunneling problems between sectors of different topology. 



Therefore, onl y a few groups have attempted to implement dynamical overlap fermions [U 

H, EE B3, G31iE El, 0, Eil, • 
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On the other hand, GW-type operators, fulfilling the GW condition in some approximation, 
although more expensive than simple Wilson-Dirac operators, have been studied in quenched cal- 
culations within the BGR collaboration for some time. There we demonstrated that at least for 
baryon masses the 0(a 2 ) corrections are quite small 23[] and that field renormalization constants 



behave almost like in the chirally symmetric case |24j . Motivated by these results we have started 
to implement CI fermions for a dynamical simulation on smaller lattices |25| | and are now presenting 
details and results of our simulations on larger lattices. 



First results involving dynamical CI fermions on 16 J x 32 lattices were published in [261 . |27| |. 



and results for smaller lattices can be found in 25|, |28|]. In this paper we concentrate on technical 
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aspects of the simulation and present first results for the hadron mass spectrum for three sets 
of parameters, corresponding to three different pion masses, giving an overview of the current 
project status. We start with an explanation of all the technicalities, i.e., simulation details and 
equilibration behavior, followed by the first analysis results for ground state masses of mesons and 
baryons. We finish with a discussion of the results and a summary. 



II. SETUP AND SIMULATION 
A. CI Dirac operator and action 

For the fermions we use the so-called Chirally Improved (CI) Dirac operator Dqi @,Bl3]. It obeys 
chiral symmetry only approximately, depending on the truncation in the extent of the interaction 
terms. Plugging a general ansatz into the Ginsparg- Wilson equation leads to a set of algebraic 
equations for the coefficients, which can be solved to obtain Dq\. The paths and coefficients 
used in our simulation are given in Appendix IA 1[ Whereas in the quenched simulations the Dqi 
coefficients were adapted to the values of the gauge coupling such as to have (almost) no mass 
renormalization, we now decided to use the same Dqi parameters for all dynamical runs. This 
implies an additive mass renormalization, i.e., the "mass parameter" m,Q does not give the bare 
mass directly. We adjust the value of mo such as to get suitable PCAC masses (also called AWI- 
masses since their definition comes from the axial Ward identity). The numbers will be discussed 
in more detail below. 

It was shown in a quenched calculation using Dqi that the Liischer-Weisz gauge action 
[29| produces smoother gauge configurations than the Wilson gauge action, and thus is used in our 
simulation. For completeness we also list details of the gauge action in Appendix I A 2[ 

Another important ingredient in our simulation is smearing, since the smearing procedure results 
in better chiral properties of the operator, as can be seen from the eigenvalue spectrum of the Dirac 



operator [30(j. In earlier quenched studies with Dqi so-called HYP smearing [3l|] was applied. Since 
such a smearing procedure is not differentiable and thus not well suited for Hybrid Monte-Carlo 
simulations, we decided to use the "differentiable" stout-smearing In our simulation we have 
been using one level of stout-smearing, such that the value of the plaquette is maximized. The 
stout smearing is considered to be part of the definition of the full Dirac operator. 

More recently, other suggestions for efficient differentiable smearing methods have been pub- 



lished 33], |34|. For consistency we continued to use the stout type smearing with which we started 



our study. 



B. Run parameters 

For the simulation presented here we use lattices of size 16 3 x 32 at three different values of 
the gauge coupling (3\ and the bare mass parameter mo, all of which can be found in Tab. [H The 
physical volume is always ~ 2.4 fm. The pion mass ranges from approximately 530 MeV down to 
320 MeV. The total number of gauge configurations produced, iV COI1 f , can also be found in Tab. [TJ 



C. Algorithm 

The algorithm we use for generating our gauge configurations is a Hybrid Monte-Carlo (HMC) 
[35| algorithm plus some additional features. HMC seems to be the most suitable algorithm for 
our goal. 
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Run 


mo 


ft 


ft 


ft 


™HB 


-^conf 


^acc 


A 


-0.050 


4.70 


-0.3941 


-0.06063 


V0.03 


100 


0.904 


B 


-0.060 


4.65 


-0.3899 


-0.05998 


V0.02 


200 


0.911 


C 


-0.077 


4.58 


-0.3841 


-0.05908 


V0.02 


200 


0.858 



TABLE I: The parameters for the different runs. The number of pseudofermions App = 2, the total length 
of the trajectories is 1 in HMC time units. In the 6-th column the parameter for the Hasenbusch mass 
preconditioning is given (see Sect. Ill Cl for more details). 



For the HMC we need a generalization of the Hamiltonian evolution for a system of classical 
mechanics in a ficticious HMC time to our system of fields U n ^- For that purpose we introduce 
traceless hermitian matrices P n „ G su(3) which act as conjugate momenta for the U n ^, with 
n, ji being the lattice site and the direction of the link, respectively. We now can define the time 
derivative of U n>tl as 

Un,ii — i Pn,fi Un,pt • (1) 

Then, a Hamiltonian H can be defined as 

H = \ E Tr ( p ?J + s z+ ^(^rv , (2) 

n,fi 

where S g denotes the gauge action and 4> is the pseudofermion field. The equation of motion for 
P is obtained via the relation H = 0, 

H = Y,^ (X/Am) + 4 + ^(Dt£>)-i0 = o , (3) 

n,fi 

which gives the evolution equation in HMC time P = f(U,U,P). Evaluating this function for, 
e.g., Wilson or staggered quarks is not complicated since such types of quarks involve only one link 
field U nt u connecting neighboring sites. In our case, however, paths up to length four, coming from 



Dqi, have to be considered. A more detailed description of the procedure can be found in [36|]. For 



the evolution in HMC time we used the reversible and area preserving leapfrog integration scheme. 



To be able to go to smaller quark masses we utilize Hasenbusch mass preconditioning 37[. The 
basic idea is to split the pseudofermion action into two (or more) parts, separating the small and 
the large eigenvalues of the Dirac matrix. In our case we always use two pseudofermions. The 



parameter ttibb, which amounts to an additional mass, is deduced from an educated guess |36l |. 
Using App pseudofermions, the mass shift is given by 



' (2 JV PF-Aj nin ) 1/7VpF , l<,<iVp F 

(4) 

» = N PF 



Here, A m j n is the assumed smallest eigenvalue of the Dirac matrix. 

For the inversion of D'D we use the standard conjugate gradient (CG) inverter. These in- 
versions take by far most of the computer time. Thus, several attempts were made to increase 
the performance of th is p art of our code. First of all we use a chronological inverter by minimal 
residue extrapolation [381 ] . taking into account 12 previous solutions. In Fig. [1] we plot the number 
of conjugate gradient iterations against the leapfrog step ?lf- What we see is a rapid decrease in 
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FIG. 1: The average number of needed conjugate gradient iterations is plotted against the first leapfrog 
steps for each run. 



the CG iteration number when more previous solutions become available. However, we find that 
a plateau is reached already at zlf = 5. The overhead caused by the 8 additional matrix vector 
multiplications is negligible, however. 

In a recent paper 39(] Diirr et al. presented a mixed precision inverter for the Dirac matrix. 
In order to ensure reversibility in the molecular dynamics (MD) evolution one should work with 
double precision accuracy. The method suggested there allows to iteratively improve the inversion 
accuracy, working partly with single precision and thus faster arithmetic. We choose a final ac- 
curacy of e = 10~ 7 . The gain in run-time per gauge configuration was, e.g., about 33% for run 
C. 



D. Autocorrelation time 

A measure for the statistical efficiency of an observable O is the integrated autocorrelation time 
ri n t, defined by 

1 ^ r(t) 

Tmt - g + f(0) ' [b) 
where the autocorrelation function T is given by 

r(t) = ((o(t )-(o))(o(t + t)-(o))) . (6) 

In practice, the sum © has to be truncated at some upper value i maX ) which we choose at that 
point where the autocorrelation data becomes noisy. We discuss several observables to be able to 
figure out the point of equilibration and the statistical independence of our measurements. 

On the l.h.s. of Fig. [2] we plot the plaquette values for the three runs. One can clearly see that 
the runs A and B, starting from "cool" quenched configurations, are equilibrated after roughly 
0(100) configurations. Run C does not show a significant equilibration process for the following 
reason. We started from a configuration B and slowly changed the parameters Pi and mo to the 
values of run C. This was the starting configuration of the new run sequence C. 

Another indicator of equilibrium behavior is the number of CG-steps in the final accept/reject 
step of the MD evolution, N[ nv . We show these numbers on the r.h.s. of Fig. [2j Here one can also 
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FIG. 2: L.h.s.: The spatially averaged plaquette against the HMC time, from top to bottom we plot runs A, 
B and C. The dashed line in the plots indicate a change in the algorithm: From that point on we changed 
to the mixed precision inverter and used Hasenbusch mass preconditioning. In run C the Hasenbusch mass 
preconditioning was used from the beginning, wc only changed to the mixed precision inverter. The full 
lines in run A and run B indicate a split in the particular run into two separate trajectories to increase 
the production of gauge configurations per (real) time. R.h.s.: The number of CG iterations, N- mv in the 
accept /reject step, notation like l.h.s. 



Run 




Tint(plaq.) 


T int(AW) 


A 


100 


3.5 


4.2 


B 


115 


2.4 


2.7 


C 


50 


3.7 


3.6 



TABLE II: Integrated autocorrelation times for the three runs. N e q U \ is the number of configurations skipped 
after the start. 



see that the runs are equilibrated after the above mentioned number of HMC updates. Based on 
these observations, in our analysis we discarded the first 100, 115, and 50 configurations for runs 
A, B, and C, respectively. 

Starting thus with the equilibrated configuration, we computed the integrated autocorrelation 
time Ti n t for the plaquette values and for -/V; nv . The resulting numbers are given in Tab. [XT] and 
are all below 5. Therefrom we decided to analyze every 5-th configuration, i.e., configurations 
separated by 5 units of HMC time. 

In Fig. [3] we show the history of the pion mass, calculated separately for each analyzed con- 
figuration. Details of the used pion interpolator are discussed in the spectroscopy section. No 
noticeable correlation can be found in the plots. 
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FIG. 3: Time histories of the pion mass determined from single configurations using the interpolator u w j5d w 
(see Sect. lIV Cj) . The horizontal dashed line indicates the mass value obtained from fits to the ww propagators 
in the range t = 4 — 15 (run A, B) or t = 5 — 15 (run C). The vertical lines in runs A and B indicate a split 
into separate sequences such as to enhance statistics by parallel runs. 



E. The change in the Hamiltonian 

Since we introduced the conjugate momenta P, we describe a microcanonical ensemble of a 
classical system with a Hamiltonian H. For exact solutions of the equations of motion (MD 
equations) the Hamiltonian would be a constant of motion and the configurations all would lie on 
a surface of constant energy. Thus, each created configuration would be accepted. However, due 
to the discretization with an MD time step 5t numerical errors are introduced and the Hamiltonian 
energy is not invariant. We denote the change as AH. Each calculated gauge configuration is then 



accepted with a probability e . The area preserving property of MD leads to an inequality [40j] , 



e (-AH) < ( e -AH^ = 1 (?) 

Due to this inequality (AH) has to be positive and this is indeed the case in our simulations 
(cf. Tab. IIIII and Fig. [5] for our values). For run B we have a quite large value of (AH), coming 
from a huge spike in AH in configuration 730 which is of the order of 0(1O 5 ) bigger than the rest. 
Also in run C we have a spike at configuration 850, being about 0(1O 3 ) bigger than the other 
values. Such spikes have already been observed in other simulations with dynamical fermions 
lil . 41]. Two possible reasons can cause such a spike. One is the instability of HMC for large step 



sizes in the MD evolution, cf. Ref. [421 ] . The other one, and this is most likely the case here, is that 
the Dirac operator can develop very small eigenvalues which lead to these spikes in the derivative 
of the action. 
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Run 


(AH) 


e -(AH) 


(e- AH ) 


A 


0.038(11) 


0.963 


0.989(11) 


B 


2.01(1.95) 


0.134 


0.986(10) 


B' 


0.055(10) 


0.947 


0.988(10) 


C 


0.089(59) 


0.915 


1.034(12) 



TABLE III: Averages of AH and their exponentials for each run. We only included the equilibrated config- 
urations in our calculations. The 3-rd row contains the data of run B without including configuration 730, 
which is responsible for the spike in AH. 
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FIG. 4: We plot AH against the HMC time starting from the point of equilibration (runs A, B and C are 
ordered from top to bottom). 



We want to conclude with a remark on the relation between AH and the acceptance rate. In 
Fig. [5] we plot the acceptance rate against the averaged AH. In our case, at least run A and run 



C are lying (within error bars) on the predicted curve [43[] 



Pace = erfc (^^j , (8) 



where erfc is the complementary error function. 
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FIG. 5: Pace vs. (AH). The black line corresponds to erfc(\/AH/2). 




FIG. 6: The smallest 150 eigenvalues superimposed for 20% of the configurations of run A. 



F. Dirac eigenvalues 



An indicator of the "GW quality" of the Dirac operator is its eigenvalue distribution in the 
complex plane. Whereas Dirac operators obeying the GW condition in its simplest form have 
eigenvalues on a unit circle centered at 1, approximate GW operators like Dqi deviate from that 
simple shape showing some scattering of the eigenvalues. Figure [6] shows the (in absolute value) 
smallest 150 eigenvalues superimposed for 20% of the configurations of run A. Obviously the fluctu- 
ation is predominantly towards values inside the unit circle and so-called exceptional configurations 
(exceptionally small eigenvalues) are suppressed. 

Figure [7] shows histograms for the smallest values of purely real A and for the minimal Re(A) for 
all three parameter sets. Both types of histograms give an indication on the permissible values of 
the smallest quark mass we may obtain for that action, lattice spacing and lattice size. Concerning 
exceptional configurations, we find a mass gap indicating that we are in a safe region of parameter 
values. 

Several observations can be made from the eigenvalue distributions. Low lying eigenvalues are 
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FIG. 7: Histograms for the smallest values of real A (left) and the smallest values of Rc(A) (right) for 
parameter sets A-C (from top to bottom). The measured AWI-mass (in lattice units) is indicated by the 
vertical, dashed line. 



depleted as expected for dynamical fermions due to the effect of the determinant in the measure. 
The boundary close to the circular shape is rather sharp towards larger values of |A — 1|. This 
allows to simulate smaller pion masses on coarse lattices. 

The distribution density towards the inner region is, for a given scale parameter, narrower than 
that for the Wilson action but not as close to the boundary as for quenched simulations with Dqi 
0. 

In the quenched simulation hypercubic smearing was used whereas for the dynamical simulation 
we apply stout smearing. This latter type of smearing has a weaker smoothing effect than the 
hypercubic type. We could have applied several subsequent stout smearing steps instead, but we 
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FIG. 8: History of the topology sector v and corresponding distribution histogram for parameter sets A-C 
(from top to bottom) 



did not want to change the effective action in the middle of our runs. Also, for the quenched 
ensembles we optimized the action parameters for each value of /3±. In the dynamical simulation 
we stayed with the same parameterization of the Dei (except for the bare "mass" parameter mo) 
in order to be able to qualitatively compare different runs. 

The number of exactly real modes u, counted according to their chirality (V'ItsIV'}) may be 
related to the topological charge via the Atiyah-Singer index theorem 44|. Although we cannot 
exclude that we miss some of the inner real modes (cf., Fig. [7]), we still get some information on the 
tunneling between topological sectors from this quantity. Figure [8] demonstrates frequent tunneling 
and consistency with a Gaussian-like shape of the distribution. 
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FIG. 9: Fits to the potential in the range 1 < r < 7 (symbols represent data points, whereas full black lines 
are fits to these points). The dashed black line in each plot indicates the distance r in lattice units where 
Eq. |n|) holds. 



III. RESULTS FOR LOW ENERGY PARAMETERS 



A. Setting the scale 



For the determination of the lattice spacing we used the Sommer parameter [45(, determined 
by the lattice potential, which was derived from Wilson loops W(r,t). For improving the signal 
the gauge configurations have been smeared with hypercubic blocking 31[ with parameter values 
a\ = 0.75, ci2 = 0.6 and 03 = 0.3. 

We have extracted the potential V(r) for each value of r from linear fits to ]n.W(r,t) in the 
range 4 < t < 7. The potential was then fitted in the range 1 < r < 7 to 



V(r)=A-\ \-ar + CAV(r) with AV(r) 



(9) 



(all quantities given in lattice units). The perturbative lattice Coulomb potential [1/v] serves as a 
correction to the continuum Coulomb potential as discussed in 46j, |47J, |48|, |49j . It has been used 
in the form corrected for hypercubic blocking 5(3, 51]. 



7T 



d 3 k cos(fc-r) -5 H Yp(fc) 
. w (2tt)3 4Ei=ism 2 (fei/2) 



(10) 



The smearing factor Shyp(&) is detailed in [50j]. The correction term allows for a perfect fit, even 
including the r = 1 value, see Fig. GO Actually, as observed by other authors, the result lies very 
close to what one gets when fitting only the continuum shape of the potential to a restricted range 
2 < r < 7. 

From the resulting potential without the correction term AV and the condition 



dV(r) 



dr 



r=ro 



we obtain the Sommer parameter in lattice units, 



The lattice spacing is thus given by a 
spacing are given in Table ITVT . 



1.65 + 5 



1.65 



^0,exp 



a a 
n),ex P /V Using 



(11) 



(12) 

0.48 fm, our values for the lattice 
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Run 


a [fin] 


aAo.cxp 


am AWI 


m AWI [MeV] 


A 


0.1507(17) 


0.3139(35) 


0.0327(3) 


42.8(4) 


B 


0.1500(12) 


0.3126(24) 


0.0259(2) 


34.1(2) 


C 


0.1440(12) 


0.3000(24) 


0.0111(2) 


15.3(3) 



TABLE IV: Lattice spacing as denned via the Sommer parameter and AWI-mass in lattice units and in 
physical units via that scale setting. 



The physical value of ro iOX p for our situation (two mass degenerate quarks) is not accessible. 
Often the scale is set by extrapolating the measured values of the lattice spacing to vanishing quark 
mass and using the extrapolated value for all mass values [52| (mass independent scheme). In the 
present state of our simulations we have only one mass value for each gauge coupling. We therefore 
rely on the mass dependent definition, which differs by 0(a) corrections. We also could use the 
nucleon mass to set the scale. In some of the mass plots shown below we therefore plot the masses 
in units of the nucleon mass. 



B. The axial Ward identity mass 

Another important observable in lattice QCD calculations is the (unrenormalized) quark mass 

from the axial Ward identity and the PCAC relation, the so-called AWI-mass (or PCAC-mass). 
Therefore, we compute the ratio 

1 c A (d t A 4 (p = 0,t)P(0)} 

mAWI = 2 ~ p W=0) t)P(0)> ' (13) 

Both interpolators 

At = (i74 75U , P = d~f 5 u, (14) 

couple to the pseudoscalar meson channel (here, time is direction 4). In relating the lattice mea- 
surements to the MS-scheme, these operators are usually defined with point-like quark sources. 
The normalization factors ca, cp relate the smeared source lattice operators to the point source 
lattice operators, 

(X^(t)P^(O)) 
X (X(")(t)PiP)(0)) ' 1 j 

where the upper index (s) or (p) indicates smeared or point sources, respectively, and X refers to 
A4 or P. The ratio c^/cp is read off from the plateau range as exhibited in Fig. [10] for the case of 
the wide sources. 

For the ratio in Eq. (|13|) we need derivatives of the correlator with respect to t. We obtain 
the numerical derivatives from local 3-point fits to the expected cosh-behavior of the correlator, 
involving values at (t — 1, t, t + 1). 

In Fig. [10] we show the AWI-mass ratio Eq. (|13p vs. t and give the corresponding numbers in 
Table IIVI The values are symmetrized with regard to T/2 and the error is estimated by single 
elimination jackknife. To obtain the final value for toawl the ratio was averaged from t = 4, . . . , 16, 
weighted according to the statistical errors. 

To relate this lattice value of toawi to an MS-value (e.g., at /i=2 GeV) we still need to extrap- 
olate to the chiral limit and compute the corresponding renormalization constants of the axial and 
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FIG. 10: L.h.s.: Ratio ca/cp for each run. R.h.s.: The AWI-mass ratios from Eq. (fT3|) (runs A, B and C 
from top to bottom). 



the pseudoscalar operators, 

m MS = m AWlZA/Zp . (16) 

The Gell-Mann-Oakes-Renner (GMOR) relation establishes (in leading order in the quark mass) 
the connection between the pion mass and the quark mass m, with F n and S denoting the 
pion decay constant and the chiral condensate, respectively, 

F>^ = -2mE. (17) 

From Fig. [11] one sees that the expected linear dependence of m 2 on ?tt,awi is nicely reproduced. 
In addition to the three fully dynamical points we also show the partially quenched values, where 
the valence quark mass is larger than the sea quark mass. These points, including the partially 
quenched ones, are all compatible with a common behavior. 



C. Pion decay constant 

The pion decay constant F n can be extracted from the correlator {A 4 A 4 ) via 

c 2 A Z A (A 4 (p = 0,t)A 4 (p = 0,0)} ^ m^e-^. (18) 

Here, we also use the normalization factor ca from Eq. (|15p in order to remove the dependence on 
the quark smearing. The axial vector lattice field operators have to be multiplied with normal- 
ization constants Z A in order to ensure correct current conservation in the chiral limit. For the 
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FIG. 11: Gcll-Mann-Oakcs-Rcnncr plot for the three runs. Full symbols represent the fully dynamical points 
whereas open symbols are data points for which m va \ > m S ea- The curves represent fits to am + bm 2 . 




0.01 0.02 0.03 0.04 
am Awi 



FIG. 12: F n in lattice units with a linear fit to the three dynamical points. These values have not been 
corrected by multiplication with Za- 



quenched results these were determined for Dq\ in 24J resulting in values close to 1 . In Fig. [12] we 
plot F n vs. the AWI-mass. 



IV. RESULTS FOR THE HADRON GROUND STATES 



A. Spectrum analysis: Variational method 



Over the last two decades lattice QCD has turned into a powerful tool for computing the mass 
spectrum of hadrons. Such a reproduction of experimental evidence from an ab-initio calculation 
is a strong test for the correctness of QCD. However, one mostly is restricted to the ground 
state masses, since excited state contributions only appear as sub-leading terms in the Euclidean 
correlators. Thus, a reliable separation of excited and ground states, but also of different excited 
states themselves, is a rather challenging enterprise. 

Nowadays several different approaches towards that goal are used in hadron spectroscopy. One 
could do a brute- force least-squares fit to a finite sum of exponentials, but this is known to give 
conclusive results only if high statistics are available. Other methods are based on Bayesian fitting 
[53L 54 . 55 . 56 . 57], subtractions [5^] or evolutionary fitting methods [o^ . 60]. Here, however, we 
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use a different state-of-the-art approach, namely the variational method 471. 161 1 which has been 
used quite extensively within the BGR collaboration (6^, 63, 64, [sB, 66, 67, 68, [6(J, 70]. For a 
recent review on results for the variational method see [711 ] . 

In the variational method a matrix built from different correlators is used. These correlators 
contain interpolators with different Dirac structure and quarks smeared with different widths. Such 
a choice allows for a better overlap of the interpolating fields with the physical states. Given a set 
of N basis interpolators Oi,i = 1, . . . , N, we compute a matrix of cross correlations, 



dj{t) = (Oi(t)Oj(0)) . 
Considering the generalized eigenvalue problem, normalized at some time slice to < t, 

C(t) • v k = X k (t, t ) • C{t ) ■ Vk , 



one can show along the lines of 6l|, [72| , that the eigenvalues behave as 

Ajfc(Mo) « e -(*-*°W [l + O ( e -(t-t0)Am fc j _ 



(19) 



(20) 



(21) 



In general, Am^ is the mass difference to the closest lying state. For a more detailed discussion of 
the error terms see [72I ] . Each of the interpolators Oj has the quantum numbers of the corresponding 
hadron channel and is projected to a certain spatial momentum, which is always zero in our case. 
For all considered hadron channels we use to = 1. 

For a sufficiently large set of basis interpolators each eigenstate decays exponentially with its 
energy according to Eq. (|2ip . The eigenstate with slowest decay (i.e., the largest eigenvalue) 
corresponds to the ground state, the second largest to the first excited state, and so on. We now 
can fit the states by stable two parameter fits of the eigenvalues in a range of i-values where the 
correlator is dominated by a single exponential. In order to identify the corresponding range for 
the fit we analyze effective masses for the eigenvalues, 



(off) 



(t + 1/2) = In 



( A fc (t) 

\A fc (t + i; 



(22) 



For sufficiently large values of t the effective masses form plateaus, which then give us the range 
for the fit. 

Another important instrument to estimate the quality of the signal are the eigenvectors of 
Eq. (|20p . acting as fingerprints for each state. The components also should show a plateau behavior 
with regard to the correlation distance where the channel is dominated by a single state. Thus, 
the fits of the eigenvalues should only be performed in a i-range where both effective masses and 
eigenvectors show a reliable plateau. 



B. Jacobi smearing of quark sources 



Hadron correlation functions are built from quark propagators D~ l acting on some quark source 
S. In order to improve the signal and to extend the operator basis we work with extended sources 
obtained by Jacobi smearing 73], [74|: A point-like source So is smeared out by acting with a 
smearing operator M, 



v 



S = MS< 







M 



(23) 
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Meson J c Number Operator 



Pseudoscalar 0^ 


1 


u, 


Ibdn 




2 


~, 


lbd w 




3 


u w "/*>d w 




4 


u, 


Itlbdn 




5 


77, 


Itlbdw 




6 




Vector 1 — ~ 


1 


a, 


"fkdn 




2 


77, 


lkd w 




3 


u w "fkd w 




4 


u, 


Ikltdn 



5 Un7kJtd w 

6 Tiwjkjtd w 



TABLE V: Meson interpolators used in this study. We use = 74, i.e., the 4-direction corresponds to the 
Euclidean time direction. The subscripts n or w denote the narrow or wide smeared quark source. 



where H is a hopping term, 



H = J2 [UM t) S s+l , + U}{x-l t) 6, : . 



(24) 



3=1 



The smearing extends only over individual time slices, i.e., t is fixed. The parameters k (hopping 
parameter) and (number of smearing steps) are tuned to get an approximately Gaussian shape 
of the quark source with a certain width. We use the values for k and N given in [(38| for the 
16 3 x 32 lattice to obtain a narrow (index n) and a wide (index w) source. 



C. Hadron interpolators 

Working with the variational method one strives for a good basis of interpolators Oi from which 
one can obtain a combination coupling strongly to the hadron of interest. These interpolators 
should simultaneously be linearly independent, as orthogonal as possible and sufficient to represent 
the physical states reasonably well. Thus, the crucial point is the design of different interpolators. 

A complete list of our meson interpolators can be found in Tab. IVl All considered interpolators 
represent isovector (/ = 1) mesons. 

Interpolators for baryons are slightly more complicated since there are three quarks involved. 
The general form of a local interpolator for the nucleon is given by 

On = e a bc Ti u a (uj T 2 d c - dl T 2 u c ) , (25) 

where a,b,c are color indices and Ti,^ are combinations of 7-matrices. In Tab. IVII the different 
possibilities are listed. 

The delta baryon has a simpler structure, since there only one Dirac structure is analyzed. Its 
interpolator has the following form, 

OA,k = tabcu a (ul Cj k u c ) , A; = 1,2, 3. (26) 



We project this to spin | and average the correlators as discussed in 691 ] . 
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r 2 


Number 


Smearing 


1 


C75 


1 


71 


ThTb) 






2 


'fi 


T)H1l 






3 


n 


HIT) 1 






4 


" \ 


711111 1 
WW) 






■5 


IV 


WTlj 






6 


w 


ww) 


it 


C7475 


13 


n 


nri) 






14 


n 


'nw) 






15 


n 


[wn) 






1G 


n( 


ww) 






17 


w 


wn) 






18 


w 


ww) 



TABLE VI: Nucleon I{J P ) = \ interpolators. The reference numbers of the interpolators are chosen 

to be consistent with earlier publications 



Number Smearing 

1 n(nn) 

2 n{nw) 

3 n(wn) 

4 n{ww) 

5 w(wn) 

6 w(ww) 



TABLE VII: Delta baryon I{J P ) = f interpolators. The reference numbers of the interpolators arc 

chosen to be consistent with earlier publications 6|| 69| . 



We introduce a short-hand notation for the different baryon interpolators. We denote them by 
si(s2S3), where Sj represents the smearing type of quark i; e.g., in n(ww) the first quark has a 
narrow smearing, the second and third have wide smearings. The interpolators for the baryons 
studied here can be found in Tables IVII and IVIIl All baryon correlators are projected to definite 
parity. 

For subsequent configurations the quark sources (and thus the hadron interpolators) are placed 
at alternating positions di = (t, x) with 

dx = ( 0, 0, 0, 0) , 

d 2 = (16, 0, 0, 0) , 

d 3 = ( 0, 8, 8, 8) , 

d 4 = (16, 8, 8, 8) , (27) 

in order to get better statistical decorrelation of the data. All hadron interpolators are projected 
to vanishing spatial momentum. 
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Meson 


Interpolator^ 


Run 


Fit range 


Mass [MeV] 


Pseudoscalar 


3 


A 


4-15 


526(7) 






B 


4-15 


469(4) 






C 


5-15 


318(5) 


Vector 


4,5,6 


A 


3-10 


922(17) 






B 


3-13 


897(13) 






C 


4-9 


810(28) 



TABLE VIII: Here we show the interpolators entering the final analysis and the best fit ranges for the 
different runs. We also give the resulting mass values using the lattice spacing given in Table HVl 



D. Effective masses and fit ranges 



Only a posteriori one can judge on the amount of independence of the interpolators used. 
Including too many interpolators in the correlation matrix increases the statistical noise in the 
diagonalization. Our aim is to get the best signal in each channel and this is obtained by having 
the best plateaus in the effective masses. We therefore, after studying the quality of results with 
different subsets of interpolators, decided on as few of them as seemed sufficient for a stable signal. 
For example, for the positive parity nucleon we only included the interpolators 4—6 and 16—18. The 
optimal selection may be different when we study higher excitations and may include differently 
smeared quark sources. 



E. The meson sector 



Since we simulate two mass-degenerate light quarks, interpolators of the form u n Vd w and u w X^d n 
are identical. Tab. [V] lists the interpolators used. Due to the two different possibilities for T we 
have for the pseudoscalar and vector particle six interpolators at hand. Only a subset of these is 
used for the final analysis. 

We restricted ourselves to fit only plateaus with three or more consecutive points. In addition 
to that we started fits only at points for which t — to > 2. Table IVThTl gives the information on the 
interpolators and the fit ranges used in the final analysis. 



1. The pseudoscalar meson 



Let us start our discussion with the particle where the best signal can be extracted, the pseu- 
doscalar meson (J = '"). For the determination of the ground state we used one interpolator 
(no. 3) and performed a cosh-fit. 

Fig. [T3l demonstrates a peculiarity of the generalized eigenvalue problem, as it was observed 



already in, e.g., Refs. [7fj, [75|. On the periodically closed lattice mesons propagate forward and 



backward in time. An interpolator which couples to a particular state at small t will also couple to 
the same, but backward running, state at high t. In the standard eigenvalue problem, depending on 
the time extent and the masses of ground state and excited state, above some value of < t\ < nt/2 
the backward running ground state will have a larger eigenvalue than the first excited state. In 
that region of t- values the second largest eigenvalue increases towards nt/2. In the generalized 
eigenvalue problem the eigenvalues are all normalized to unity at timeslice to- Thus the second 
largest eigenvalue signal is shifted upwards and the upwards increasing eigenvalue discussed may 
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FIG. 13: In the first row the eigenvalues for the pscudoscalar channel (J PC = h ) arc shown (runs A. B, 
C from left to right). In each plot we show data for two different sets of interpolators: Circles represent the 
ground state using interpolator 3, squares and diamonds show the ground state (GS) and the first excited 
state (IE) using the interpolator set 1, 3, 4, 6, respectively (numbers according to Tab. lV)) . In the second row 
the absolute value of the corresponding effective masses (in lattice units) of the ground states are plotted as 
a function of t. The horizontal line indicates the fit range and mass value obtained by the fit of the ground 
state eigenvalue of interpolator 3 over the specified range. 



now even becomes at some value ti larger than the eigenvalue of the ground state. This behavior 
is observed in Fig. [13] where we plot the first two eigenvalues of the pseudoscalar state resulting 
from the generalized eigenvalue problem analysis. For our choice of to = 1 the (in time) backward 
running ground state becomes the second largest eigenvalue near t\ = 6 and the largest eigenvalue 
near t<i = 13. Near the crossing this leads to a misidentification (the real ground state signal 
becomes the second largest eigenvalue) which explains the bump in the effective mass |am e g|. 

These properties can be seen very nicely in the 2-dimensional model of [7g]. This behavior 
of the eigenvalues is a fundamental feature of the variational method; the signals of ground and 
excited states are disentangled up to that point in time where these signals are crossing with the 
lightest backward running state. The larger the difference in the ground and excited states, the 
earlier this crossing takes place. 

For simplicity, and since the results for the corresponding plateau regions agree within errors, 
we choose the single correlator value where we find the longest plateau. In Fig. [13] we compare the 
effective masses of the ground states for two different choices of interpolators. One can clearly see 
that the two sets of effective masses can be fitted reliably in an appropriate region. 
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FIG. 14: In the first row the normalized eigenvector components of the ground state 



4, 5, 6, of the 



vector meson (J PC = arc plotted against the time distance t. From left to right we show runs A, B 

and C. In the second row the effective masses of the vector meson (in lattice units) are plotted as a function 
of t. The solid black lines indicate our fit range of the fit to the corresponding leading eigenvalue and the 
upper and lower bound of the extracted value. Also here, from left to right we show runs A, B and C. 



2. The vector meson 



In the case of the vector meson we can take the eigenvector components as a tool to determine 
fit ranges (cf., Table Em]). The interpolators we included in the correlation matrix are no. 4, 5, 6. 
In Fig. Q3] we plot the eigenvector components and effective mass of the ground state. One can see 
from the plots that the quality of the data is sufficient to make a fit, but it is not as good as for 
the pion. 

On the l.h.s. of Fig.[15]we plot the fitted mass against m^. The scale is set by the lattice spacing 
of Table IIVI Within error bars, all three runs agree nicely with each other and run C extrapolates 
close to the experimental value. 

As discussed in Sect. IIIII we set the scale by assuming as Sommer parameter value of 0.48 fm 
for all three runs. In the r.h.s. of Fig. [15] we use the scale of the nucleon mass instead. Here the 
data of run B seem to be somewhat higher than runs A and C. 

We emphasize that the physical p is a resonance and that multiple lattice volumes would be 
needed for a thorough analysis. 



F. The baryon sector 

In this presentation we restrict ourselves to baryons with positive parity. A more detailed 
analysis, also including excited states, is in progress. The definitions for the nucleon and delta 
baryons were given earlier in Eqs. (|25p and (126ft . Details of the interpolators used can be found in 
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FIG. 15: L.h.s.: The vector meson mass my is plotted against the pseudoscalar mass R.h.s.: An APE 
plot (scaled by the nucleon mass m^)- Both plots show the fully dynamical data (filled symbols) and the 
partially quenched dynamical data (open symbols) for runs A, B, C. The physical point is marked with a 
black cross. 
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FIG. 16: The (positive parity) nucleon mass m^r is plotted against the pseudoscalar mass for the fully 
dynamical data (filled symbols) and for the partially quenched dynamical data (open symbols) for runs A, 
B, C. The experimental value is marked with a black cross. 



Tables ED and Em 



1. The nucleon 



For the diagonalization process we used the interpolators 4, 5, 6, 16, 17, 18 according to Tab. IVIl 
The results are shown in Fig. [TBI All of the data sets extrapolate towards the physical point. 



2. The delta resonance 



For the A resonance we have a set of 6 different interpolators at hand (see Tab. IVIip and in 
principle we can allow for 2 6 — 1 = 63 combinations. All these combinations give rise to reasonable 
fit results. In the end we used the combination 1, 2, 4, 5, 6, see Fig. [TTJ However, a naive (linear in 
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FIG. 17: L.h.s.: The (positive parity) delta baryon mass toa is plotted against the pscudoscalar mass m^. 
R.h.s.: An APE plot (scaled by the nucleon mass) for the dynamical runs A, B, C. Both plots show the fully 
dynamical data (filled symbols) and the partially quenched dynamical data (open symbols) for runs A, B, 
C. The physical point is marked with a black cross. 



Baryon 


Intcrpolator(s) 


Run 


Fit range 


Mass [McV] 


Nucleon (pos. parity) 


4,5,6,16,17,18 


A 


3-11 


1311(22) 






B 


4-11 


1215(18) 






C 


3-8 


1108(23) 


Delta (pos. parity) 


1,2,4,5,6 


A 


3-6 


1528(22) 






B 


3-6 


1498(15) 






C 


3-6 


1443(23) 



TABLE IX: Interpolators and fit ranges used for the baryon ground states. The mass values are obtained 
using the lattice spacing given in Table |lVl 



m^) extrapolation overestimates the physical value by about 10% - 15%. On the r.h.s. of Fig. [T7] 
we show an APE plot, scaled by the nucleon mass. One may argue that in this plot some finite 
size artefacts cancel such that the extrapolation to the physical point is improved. 



V. SUMMARY AND CONCLUSIONS 

In this paper we presented first results from dynamical simulations with CI fermions on lattices 
of size 16 3 x 32 with spatial extent of 2.4 fm. After detailing the technical aspects of our simulation 
we showed that so-called exceptional configurations are suppressed in simulations with CI fermions. 
This enables us to simulate at pion masses of roughly 320 MeV on rather coarse lattices. We observe 
frequent tunneling between topological sectors and reasonably small autocorrelation times. 

As a first physical application we presented results for the pion decay constant F n and for the 
ground state masses of selected mesons and baryons. While scale setting remains an issue with 
dynamical simulations, the results from all three runs are consistent and naive extrapolations of 
our data are also consistent with experiment. Further simulations at different lattice spacings and 
in larger volumes will be needed in order to control the effects of the lattice discretization and to 
estimate the finite volume corrections, thereby making closer contact with experimental results. 

We are currently improving the basis for the variational method and investigating the effects of 
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quark and link smearing on the quality of excited state signals, thus providing a systematic study 
of excited meson and baryon states for a larger set of quantum numbers. 
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APPENDIX A: CI OPERATOR AND LUSCHER-WEISZ GAUGE ACTION 

1. The CI operator 



Throughout the dynamical simulations we used the CI Dirac operator introduced in 0, 77]. 
The coefficients multiply terms of the action according to the definition in 
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D = m t + Da , D C i(n, m) = J2 c nl( U ) r * > ( A1 ) 

i=i 

where the sum runs over all 16 elements of the Clifford algebra. To each element we assign a 
coefficient Cnm, consisting of sums of path ordered products of the link variables U which connect 
the lattice sites n and m. Plugging this ansatz into the Ginsparg- Wilson equation leads to a set 
of algebraic equations, which can be solved to obtain Dq\. Additional restrictions come from the 
lattice symmetries and 75-hermiticity. The solution can in principle be exact if one allows for an 
infinite number of terms. For practical reasons the number of terms is finite and thus the solution 
is a truncated series solution of the Ginsparg- Wilson relation. In our simulation paths up to length 
four are used, given in Table [X]. 



2. The Liischer-Weisz gauge action 

The Liischer-Weisz gauge action (2f| is given by 

s z = -foJ2l Retr u & - & 1 Retr u ** - & \ Retr Uth ' ( A2 ) 

pi re tb 

where U p \ is the usual Wilson plaquette, U rc is a planar (2 x l)-plaquette and U t b is a closed loop of 
length 6 along the edges of a 3-cube ("twisted bent"). Here, /3i is the independent gauge coupling 
and the two other couplings are determined from tadpole-improved perturbation theory 78J]. With 

f\ \V4 1 

u = - ReTr (f7 Dl ) , a = logi^ , (A3) 

y 3 \ P i/ j , 3.06839 y J 

we get for ^2,(^3 the following expressions, 

fa = — (1 + 0.4805a) , /? 3 = % 0.03325a . (A4) 
20uq Uq 
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Coeff. number Name Value Path shape 7 Multiplicity 



1 


Sl 


"1 1 O 1 p" A A A r" A 

1.481599252 


r 1 

[ ] 
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1 
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S2 


— 0. 05218251439 


r -i 


1 


8 


3 


S3 


—0.01473643847 
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48 


5 


S5 


A AAA 1 0/^1 AO I A 1 

— 0.002186103421 


r ■ - 7 i 


1 


192 


G 


S6 


A A A A 1 O O A A A /"> A /"> 

0.002133989696 


r ■ -i 


1 


96 


8 


S8 


A AAOAA^TAAI AA1 

—0.003997001821 




1 
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10 


SlO 


A A A A , j A p" 1 /'iWOno P- 

—0.0004951673735 


r * ■ 7 71 

[i,3,k,l\ 


1 


384 
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Sll 


A A A A A A O r* A A ^7 A A 

—0.0009836500799 


r * -7i 

l«)j,-«,«J 


1 


384 


13 


Sl3 


0.007529838581 


r ■ ■ - i 

[h3>-h-3\ 


1 


48 


14 


('1 


0.1972229309 


r 'i 


7« 


8 


15 


''2 


a a a o o r o i r^rrr 

0.008252157565 




Ti 


96 


1 1 


Vi 


n aaci 1 oncroi /I 

0.005113056314 


r ■ ■ / 1 


7i 


OO A 

384 


18 


I'F, 


n AA1 ^orrAA/iot: 

0.001736609425 


r ■ ■ / 1 


7i 


192 


32 


h 


-0.08792744664 


M 


7*7^ 


48 


33 


to 


-0.002553055577 




7i7j 


384 


34 


h 


0.002093792069 




7i7i 


192 


36 


t-> 


-0.005567377075 


[ij, ~i] 


7i7j 


48 


4G 


1 15 


-0.003427310798 


[j, i, -3, -i] 


7i7i 


48 


51 


Pi 


-0.008184103136 


[i,3,k,l] 


75 


384 



TABLE X: Coefficients for the CI fermion action used in this simulation. The path shapes are given 
symbolically, e.g., stands for a path in i-direction and then in j-direction (i ^ j). The 7-matriccs (5-th 
column) are also permuted as described in more detail in [Io| ]. 

By uq in Eq. (|A3|) we denote the assumed plaquette, ReTr (U p i), thus the coefficients have to be 
calculated self-consistentiy. 
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